Chapter 10 Changes after FMT compared to acclimation
10.1 What is the trend of the microbiota in each type after FMT?
10.1.1 Alpha diversity
label_map <- c(
"Control" = "Cold-control",
"Treatment" = "Cold-intervention",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity"
)
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point %in% c("FMT1","Acclimation", "FMT2")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2") Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(5.7695) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
757.7 783.7 -367.8 735.7 68
Scaled residuals:
Min 1Q Median 3Q Max
-2.01194 -0.49654 0.06198 0.55787 1.90952
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.02586 0.1608
Number of obs: 79, groups: individual, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9089 0.1696 23.046 <2e-16 ***
typeHot_control 0.5568 0.2283 2.439 0.0147 *
typeTreatment -0.1233 0.2297 -0.537 0.5916
time_pointFMT1 -0.3140 0.2167 -1.449 0.1473
time_pointFMT2 0.1487 0.2173 0.684 0.4937
typeHot_control:time_pointFMT1 -0.1155 0.2983 -0.387 0.6987
typeTreatment:time_pointFMT1 0.4309 0.3068 1.404 0.1602
typeHot_control:time_pointFMT2 -0.3876 0.2988 -1.297 0.1945
typeTreatment:time_pointFMT2 0.4093 0.2990 1.369 0.1710
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typHt_cntrl -0.741
typeTretmnt -0.721 0.535
tim_pntFMT1 -0.677 0.503 0.497
tim_pntFMT2 -0.702 0.520 0.507 0.532
typH_:_FMT1 0.496 -0.671 -0.362 -0.727 -0.389
typTr:_FMT1 0.482 -0.357 -0.664 -0.707 -0.378 0.515
typH_:_FMT2 0.515 -0.684 -0.370 -0.388 -0.730 0.516 0.276
typTr:_FMT2 0.498 -0.370 -0.686 -0.384 -0.719 0.280 0.511 0.524 $emmeans
type emmean SE df asymp.LCL asymp.UCL
Control 3.85 0.1050 Inf 3.65 4.06
Hot_control 4.24 0.0996 Inf 4.05 4.44
Treatment 4.01 0.1030 Inf 3.81 4.21
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Control - Hot_control -0.389 0.144 Inf -2.711 0.0184
Control - Treatment -0.157 0.145 Inf -1.083 0.5246
Hot_control - Treatment 0.232 0.143 Inf 1.627 0.2342
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates $emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 4.05 0.0938 Inf 3.87 4.24
FMT1 3.84 0.0949 Inf 3.66 4.03
FMT2 4.21 0.0895 Inf 4.03 4.38
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - FMT1 0.209 0.123 Inf 1.700 0.2050
Acclimation - FMT2 -0.156 0.121 Inf -1.286 0.4028
FMT1 - FMT2 -0.365 0.122 Inf -2.997 0.0077
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates Neutral
Modelq1n <- nlme::lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Modelq1n)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
563.2664 587.9998 -270.6332
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.551076 9.160157
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 21.343055 3.603126 46 5.923483 0.0000
typeHot_control 23.258307 4.960549 24 4.688656 0.0001
typeTreatment -3.944072 4.960549 24 -0.795088 0.4344
time_pointFMT1 -3.983722 4.472618 46 -0.890691 0.3777
time_pointFMT2 2.374489 4.472618 46 0.530895 0.5980
typeHot_control:time_pointFMT1 -14.877819 6.216964 46 -2.393101 0.0208
typeTreatment:time_pointFMT1 11.253000 6.325237 46 1.779064 0.0818
typeHot_control:time_pointFMT2 -14.775497 6.216964 46 -2.376642 0.0217
typeTreatment:time_pointFMT2 20.108022 6.216964 46 3.234380 0.0023
Correlation:
(Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typeHot_control -0.726
typeTreatment -0.726 0.528
time_pointFMT1 -0.663 0.481 0.481
time_pointFMT2 -0.663 0.481 0.481 0.534
typeHot_control:time_pointFMT1 0.477 -0.649 -0.346 -0.719 -0.384
typeTreatment:time_pointFMT1 0.469 -0.340 -0.638 -0.707 -0.378 0.509
typeHot_control:time_pointFMT2 0.477 -0.649 -0.346 -0.384 -0.719 0.518 0.272
typeTreatment:time_pointFMT2 0.477 -0.346 -0.649 -0.384 -0.719 0.276 0.509 0.518
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.1324340 -0.7088533 0.0456110 0.6569569 1.9681295
Number of Observations: 79
Number of Groups: 27 $emmeans
type emmean SE df lower.CL upper.CL
Control 20.8 2.36 26 16.0 25.7
Hot_control 34.2 2.33 24 29.4 39.0
Treatment 27.3 2.36 24 22.4 32.2
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Control - Hot_control -13.37 3.31 24 -4.038 0.0013
Control - Treatment -6.51 3.33 24 -1.952 0.1461
Hot_control - Treatment 6.86 3.31 24 2.073 0.1170
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 27.8 2.01 24 23.6 31.9
FMT1 22.6 2.01 24 18.4 26.7
FMT2 31.9 1.97 24 27.9 36.0
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 5.19 2.55 46 2.034 0.1156
Acclimation - FMT2 -4.15 2.52 46 -1.646 0.2372
FMT1 - FMT2 -9.34 2.52 46 -3.703 0.0016
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates Phylogenetic
Model_phylo <- nlme::lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
269.7976 294.5311 -123.8988
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.4986226 1.145917
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 5.256251 0.4407611 46 11.925398 0.0000
typeHot_control 1.266332 0.6064636 24 2.088060 0.0476
typeTreatment 0.283689 0.6064636 24 0.467776 0.6442
time_pointFMT1 -0.837008 0.5590601 46 -1.497169 0.1412
time_pointFMT2 -0.481515 0.5590601 46 -0.861293 0.3935
typeHot_control:time_pointFMT1 -1.410952 0.7774021 46 -1.814958 0.0761
typeTreatment:time_pointFMT1 -0.545877 0.7906304 46 -0.690432 0.4934
typeHot_control:time_pointFMT2 -0.585100 0.7774021 46 -0.752635 0.4555
typeTreatment:time_pointFMT2 0.056099 0.7774021 46 0.072162 0.9428
Correlation:
(Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typeHot_control -0.727
typeTreatment -0.727 0.528
time_pointFMT1 -0.676 0.492 0.492
time_pointFMT2 -0.676 0.492 0.492 0.533
typeHot_control:time_pointFMT1 0.486 -0.663 -0.353 -0.719 -0.383
typeTreatment:time_pointFMT1 0.478 -0.348 -0.652 -0.707 -0.377 0.509
typeHot_control:time_pointFMT2 0.486 -0.663 -0.353 -0.383 -0.719 0.517 0.271
typeTreatment:time_pointFMT2 0.486 -0.353 -0.663 -0.383 -0.719 0.276 0.509 0.517
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.81722823 -0.52018811 -0.06481484 0.55835817 2.23628404
Number of Observations: 79
Number of Groups: 27 $emmeans
type emmean SE df lower.CL upper.CL
Control 4.82 0.280 26 4.24 5.39
Hot_control 5.42 0.276 24 4.85 5.99
Treatment 4.94 0.280 24 4.36 5.52
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Control - Hot_control -0.601 0.393 24 -1.527 0.2963
Control - Treatment -0.120 0.396 24 -0.304 0.9505
Hot_control - Treatment 0.481 0.393 24 1.221 0.4525
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 5.77 0.245 24 5.27 6.28
FMT1 4.28 0.245 24 3.78 4.79
FMT2 5.12 0.241 24 4.62 5.61
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 1.489 0.319 46 4.666 0.0001
Acclimation - FMT2 0.658 0.316 46 2.085 0.1042
FMT1 - FMT2 -0.831 0.316 46 -2.635 0.0302
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates 10.1.1.1 Beta diversity
Number of samples used
samples_to_keep_post7 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2") %>%
select(Tube_code) %>%
pull()
subset_meta_post7 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2")
subset_meta_post7$time_point<-as.factor(subset_meta_post7$time_point)
subset_meta_post7$type<-as.factor(subset_meta_post7$type)
length(samples_to_keep_post7)[1] 79Richness
richness_post7 <- as.matrix(beta_q0n$S)
richness_post7 <- as.dist(richness_post7[rownames(richness_post7) %in% samples_to_keep_post7,
colnames(richness_post7) %in% samples_to_keep_post7])
betadisper(richness_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.06454 0.032268 4.9753 999 0.01 **
Residuals 76 0.49290 0.006485
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.0160000 0.697
Hot_control 0.0119513 0.002
Treatment 0.6817089 0.0035106 adonis2(richness_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(richness_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(richness_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 2 | 3.120806 | 0.11961754 | 5.546155 | 0.001 |
| time_point | 2 | 1.703776 | 0.06530413 | 3.027874 | 0.001 |
| type:time_point | 4 | 1.570884 | 0.06021051 | 1.395852 | 0.001 |
| Residual | 70 | 19.694403 | 0.75486782 | NA | NA |
| Total | 78 | 26.089870 | 1.00000000 | NA | NA |
Neutral
neutral_post7 <- as.matrix(beta_q1n$S)
neutral_post7 <- as.dist(neutral_post7[rownames(neutral_post7) %in% samples_to_keep_post7,
colnames(neutral_post7) %in% samples_to_keep_post7])
betadisper(neutral_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.06738 0.033689 3.8309 999 0.031 *
Residuals 76 0.66833 0.008794
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.2150000 0.186
Hot_control 0.1985046 0.009
Treatment 0.1611083 0.0069718 adonis2(neutral_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(neutral_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(neutral_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 2 | 3.705559 | 0.14953527 | 7.609355 | 0.001 |
| time_point | 2 | 2.217656 | 0.08949199 | 4.553951 | 0.001 |
| type:time_point | 4 | 1.813192 | 0.07317011 | 1.861693 | 0.001 |
| Residual | 70 | 17.044094 | 0.68780262 | NA | NA |
| Total | 78 | 24.780501 | 1.00000000 | NA | NA |
Phylogenetic
phylo_post7 <- as.matrix(beta_q1p$S)
phylo_post7 <- as.dist(phylo_post7[rownames(phylo_post7) %in% samples_to_keep_post7,
colnames(phylo_post7) %in% samples_to_keep_post7])
betadisper(phylo_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.05438 0.027192 2.1499 999 0.116
Residuals 76 0.96123 0.012648
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.704000 0.164
Hot_control 0.697512 0.083
Treatment 0.155309 0.083221 adonis2(phylo_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(phylo_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(phylo_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 2 | 0.2941494 | 0.07500331 | 3.592254 | 0.001 |
| time_point | 2 | 0.5528857 | 0.14097688 | 6.752032 | 0.001 |
| type:time_point | 4 | 0.2088312 | 0.05324856 | 1.275159 | 0.147 |
| Residual | 70 | 2.8659522 | 0.73077126 | NA | NA |
| Total | 78 | 3.9218184 | 1.00000000 | NA | NA |
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_richness_nmds_post7 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_neutral_nmds_post7 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_phylogenetic_nmds_post7 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_nmds_post7, beta_neutral_nmds_post7, beta_phylogenetic_nmds_post7, ncol=3, nrow=1, common.legend = TRUE, legend="right")
10.2 What is the effect of FMT and temperature on the GM after 1 week compared to acclimation?
10.2.1 CI vs CC
10.2.1.1 Alpha diversity
label_map <- c(
"Control" = "Cold-control",
"Treatment" = "Cold-intervention",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity"
)
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(type %in% c("Control","Treatment") & time_point %in% c("FMT1","Acclimation")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control") Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(2.5327) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
329.2 338.4 -158.6 317.2 28
Scaled residuals:
Min 1Q Median 3Q Max
-1.4046 -0.7191 0.1331 0.6747 1.4330
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 2.062e-11 4.541e-06
Number of obs: 34, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9512 0.2275 17.368 <2e-16 ***
typeTreatment -0.1372 0.3132 -0.438 0.661
time_pointFMT1 -0.3225 0.3140 -1.027 0.304
typeTreatment:time_pointFMT1 0.4500 0.4435 1.015 0.310
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typTrt t_FMT1
typeTretmnt -0.726
tim_pntFMT1 -0.725 0.526
typTr:_FMT1 0.513 -0.706 -0.708
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Modelq1n <- lme(neutral ~ type*time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Modelq1n)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
244.9928 253.4 -116.4964
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 6.116057 8.521967
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 21.027978 3.684717 16 5.706810 0.0000
typeTreatment -3.628995 5.079636 16 -0.714420 0.4853
time_pointFMT1 -3.668645 4.182131 14 -0.877219 0.3952
typeTreatment:time_pointFMT1 10.917012 5.914427 14 1.845828 0.0862
Correlation:
(Intr) typTrt t_FMT1
typeTreatment -0.725
time_pointFMT1 -0.611 0.443
typeTreatment:time_pointFMT1 0.432 -0.582 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.79785802 -0.67295800 0.08387402 0.71606547 1.30935708
Number of Observations: 34
Number of Groups: 18 Phylogenetic
Model_phylo <- lme(phylogenetic ~ type*time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
125.282 133.6892 -56.64101
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 7.854299e-05 1.386163
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 5.227341 0.4900827 16 10.666243 0.0000
typeTreatment 0.312600 0.6735542 16 0.464105 0.6488
time_pointFMT1 -0.808097 0.6735542 14 -1.199751 0.2501
typeTreatment:time_pointFMT1 -0.541741 0.9525495 14 -0.568728 0.5786
Correlation:
(Intr) typTrt t_FMT1
typeTreatment -0.728
time_pointFMT1 -0.728 0.529
typeTreatment:time_pointFMT1 0.514 -0.707 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.630942795 -0.368275262 0.004976243 0.364227875 2.050663052
Number of Observations: 34
Number of Groups: 18 10.2.1.2 Beta diversity
Number of samples used
samples_to_keep_post3 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control") %>%
select(Tube_code) %>%
pull()
subset_meta_post3 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
subset_meta_post3$time_point<-as.factor(subset_meta_post3$time_point)
subset_meta_post3$type<-as.factor(subset_meta_post3$type)
length(samples_to_keep_post3)[1] 34Richness
richness_post3 <- as.matrix(beta_q0n$S)
richness_post3 <- as.dist(richness_post3[rownames(richness_post3) %in% samples_to_keep_post3,
colnames(richness_post3) %in% samples_to_keep_post3])
betadisper(richness_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.002424 0.0024241 0.4717 999 0.525
Residuals 32 0.164445 0.0051389
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.515
Treatment 0.49714 adonis2(richness_post3 ~ type*time_point,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 0.4269704 | 0.03612811 | 1.312996 | 0.001 |
| time_point | 1 | 1.1427754 | 0.09669595 | 3.514201 | 0.001 |
| type:time_point | 1 | 0.4928541 | 0.04170286 | 1.515598 | 0.070 |
| Residual | 30 | 9.7556341 | 0.82547308 | NA | NA |
| Total | 33 | 11.8182341 | 1.00000000 | NA | NA |
Neutral
neutral_post3 <- as.matrix(beta_q1n$S)
neutral_post3 <- as.dist(neutral_post3[rownames(neutral_post3) %in% samples_to_keep_post3,
colnames(neutral_post3) %in% samples_to_keep_post3])
betadisper(neutral_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.019686 0.0196859 2.8068 999 0.11
Residuals 32 0.224435 0.0070136
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.103
Treatment 0.10361 adonis2(neutral_post3 ~ type*time_point,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 0.4361022 | 0.04050485 | 1.595030 | 0.001 |
| time_point | 1 | 1.6965854 | 0.15757761 | 6.205206 | 0.001 |
| type:time_point | 1 | 0.4315814 | 0.04008496 | 1.578495 | 0.125 |
| Residual | 30 | 8.2023961 | 0.76183257 | NA | NA |
| Total | 33 | 10.7666650 | 1.00000000 | NA | NA |
Phylogenetic
phylo_post3 <- as.matrix(beta_q1p$S)
phylo_post3 <- as.dist(phylo_post3[rownames(phylo_post3) %in% samples_to_keep_post3,
colnames(phylo_post3) %in% samples_to_keep_post3])
betadisper(phylo_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04697 0.046974 2.8076 999 0.111
Residuals 32 0.53540 0.016731
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.11
Treatment 0.10357 adonis2(phylo_post3 ~ type*time_point,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 0.03710337 | 0.01720754 | 0.6366663 | 0.004 |
| time_point | 1 | 0.32272550 | 0.14967132 | 5.5377295 | 0.003 |
| type:time_point | 1 | 0.04807164 | 0.02229432 | 0.8248736 | 0.447 |
| Residual | 30 | 1.74832756 | 0.81082682 | NA | NA |
| Total | 33 | 2.15622808 | 1.00000000 | NA | NA |
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post3, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post3$time_pointFMT1" = "FMT1",
"subset_meta_post3$typeTreatment" = "Cold-intervention",
"subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")
beta_richness_nmds_post3 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post3, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post3$time_pointFMT1" = "FMT1",
"subset_meta_post3$typeTreatment" = "Cold-intervention",
"subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")
beta_neutral_nmds_post3 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post3, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post3$time_pointFMT1" = "FMT1",
"subset_meta_post3$typeTreatment" = "Cold-intervention",
"subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")
beta_phylogenetic_nmds_post3 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_nmds_post3, beta_neutral_nmds_post3, beta_phylogenetic_nmds_post3, ncol=3, nrow=1, common.legend = TRUE, legend="right")
10.2.2 Functional differences
CC from acclimation to FMT1
significant_element<-element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))element_gift_sig <- element_gift %>%
select(Tube_code, all_of(intersect(
significant_element$trait,
colnames(element_gift)
))) %>%
left_join(., sample_metadata[c(1, 7, 10)], by = join_by(Tube_code == Tube_code)) %>%
filter(time_point %in% c("FMT1","Acclimation"))%>%
filter(type=="Control")
difference_table <- element_gift_sig %>%
select(-Tube_code, -type) %>%
group_by(time_point) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-FMT1)%>%
mutate(group_color = ifelse(Difference <0,"FMT1", "Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c("#4477AA",'#008080')) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
legend.title = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Function") +
ylab("Mean difference")
CI from acclimation to FMT1
significant_element<-element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Treatment") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))element_gift_sig <- element_gift %>%
select(Tube_code, all_of(intersect(
significant_element$trait,
colnames(element_gift)
))) %>%
left_join(., sample_metadata[c(1, 7, 10)], by = join_by(Tube_code == Tube_code)) %>%
filter(time_point %in% c("FMT1","Acclimation"))%>%
filter(type=="Treatment")
difference_table <- element_gift_sig %>%
select(-Tube_code, -type) %>%
group_by(time_point) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-FMT1)%>%
mutate(group_color = ifelse(Difference <0,"FMT1", "Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c("#76b183",'#008080')) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
legend.title = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Function") +
ylab("Mean difference")
10.2.3 CI vs WC
10.2.3.1 Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=8))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control") Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(5.0429) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
338.8 348.1 -163.4 326.8 29
Scaled residuals:
Min 1Q Median 3Q Max
-1.85185 -0.65278 0.06667 0.54823 1.97668
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 1.315e-11 3.627e-06
Number of obs: 35, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.4697 0.1527 29.279 < 2e-16 ***
typeTreatment -0.6557 0.2186 -2.999 0.00271 **
time_pointFMT1 -0.4208 0.2174 -1.936 0.05292 .
typeTreatment:time_pointFMT1 0.5484 0.3146 1.743 0.08131 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typTrt t_FMT1
typeTretmnt -0.698
tim_pntFMT1 -0.702 0.490
typTr:_FMT1 0.485 -0.695 -0.691
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Modelq1n <- lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta,)
summary(Modelq1n)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
250.0982 258.7021 -119.0491
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 5.122922 8.546652
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 44.60136 3.321472 16 13.428192 0e+00
typeTreatment -27.20238 4.697271 16 -5.791103 0e+00
time_pointFMT1 -18.86154 4.028930 15 -4.681526 3e-04
typeTreatment:time_pointFMT1 26.15805 5.809237 15 4.502838 4e-04
Correlation:
(Intr) typTrt t_FMT1
typeTreatment -0.707
time_pointFMT1 -0.606 0.429
typeTreatment:time_pointFMT1 0.421 -0.595 -0.694
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.903393212 -0.678106330 0.002306491 0.652929042 1.401732098
Number of Observations: 35
Number of Groups: 18 Phylogenetic
Model_phylo <- lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
124.0635 132.6674 -56.03174
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.0002455305 1.282332
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 6.522583 0.4274440 16 15.259503 0.0000
typeTreatment -0.982643 0.6044971 16 -1.625554 0.1236
time_pointFMT1 -2.247960 0.6044971 15 -3.718727 0.0021
typeTreatment:time_pointFMT1 0.898121 0.8681429 15 1.034531 0.3173
Correlation:
(Intr) typTrt t_FMT1
typeTreatment -0.707
time_pointFMT1 -0.707 0.500
typeTreatment:time_pointFMT1 0.492 -0.696 -0.696
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.40177643 -0.54641173 -0.06022755 0.42683546 2.21670618
Number of Observations: 35
Number of Groups: 18 10.2.3.2 Beta diversity
Number of samples used
samples_to_keep_post4 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
select(Tube_code) %>%
pull()
subset_meta_post4 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control")
subset_meta_post4$time_point<-as.factor(subset_meta_post4$time_point)
subset_meta_post4$type<-as.factor(subset_meta_post4$type)
length(samples_to_keep_post4)[1] 35Richness
richness_post4 <- as.matrix(beta_q0n$S)
richness_post4 <- as.dist(richness_post4[rownames(richness_post4) %in% samples_to_keep_post4,
colnames(richness_post4) %in% samples_to_keep_post4])
betadisper(richness_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.076276 0.076276 19.518 999 0.001 ***
Residuals 33 0.128967 0.003908
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.001
Treatment 0.00010134 adonis2(richness_post4 ~ type*time_point,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 1.1360754 | 0.09995918 | 4.050611 | 0.001 |
| time_point | 1 | 0.9251516 | 0.08140076 | 3.298575 | 0.001 |
| type:time_point | 1 | 0.6095926 | 0.05363586 | 2.173467 | 0.004 |
| Residual | 31 | 8.6945735 | 0.76500420 | NA | NA |
| Total | 34 | 11.3653932 | 1.00000000 | NA | NA |
Neutral
neutral_post4 <- as.matrix(beta_q1n$S)
neutral_post4 <- as.dist(neutral_post4[rownames(neutral_post4) %in% samples_to_keep_post4,
colnames(neutral_post4) %in% samples_to_keep_post4])
betadisper(neutral_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.087335 0.087335 15.924 999 0.001 ***
Residuals 33 0.180987 0.005484
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.003
Treatment 0.00034565 adonis2(neutral_post4 ~ type*time_point,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 1.3325426 | 0.12145986 | 5.413823 | 0.001 |
| time_point | 1 | 1.2692344 | 0.11568938 | 5.156616 | 0.001 |
| type:time_point | 1 | 0.7390264 | 0.06736148 | 3.002499 | 0.006 |
| Residual | 31 | 7.6302501 | 0.69548928 | NA | NA |
| Total | 34 | 10.9710536 | 1.00000000 | NA | NA |
Phylogenetic
phylo_post4 <- as.matrix(beta_q1p$S)
phylo_post4 <- as.dist(phylo_post4[rownames(phylo_post4) %in% samples_to_keep_post4,
colnames(phylo_post4) %in% samples_to_keep_post4])
betadisper(phylo_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.07659 0.076588 5.0159 999 0.026 *
Residuals 33 0.50388 0.015269
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.027
Treatment 0.031971 adonis2(phylo_post4 ~ type*time_point,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 0.20707792 | 0.09373244 | 3.875413 | 0.011 |
| time_point | 1 | 0.28096358 | 0.12717629 | 5.258166 | 0.001 |
| type:time_point | 1 | 0.06475669 | 0.02931169 | 1.211906 | 0.289 |
| Residual | 31 | 1.65644672 | 0.74977958 | NA | NA |
| Total | 34 | 2.20924491 | 1.00000000 | NA | NA |
RDA
#Richness
cca_ord <- capscale(formula = richness_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post4, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post4$time_pointFMT1" = "FMT1",
"subset_meta_post4$typeTreatment" = "Cold-intervention",
"subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")
beta_richness_nmds_post4 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post4, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post4$time_pointFMT1" = "FMT1",
"subset_meta_post4$typeTreatment" = "Cold-intervention",
"subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")
beta_neutral_nmds_post4 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post4, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post4$time_pointFMT1" = "FMT1",
"subset_meta_post4$typeTreatment" = "Cold-intervention",
"subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")
beta_phylogenetic_nmds_post4 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_nmds_post4, beta_neutral_nmds_post4, beta_phylogenetic_nmds_post4, ncol=3, nrow=1, common.legend = TRUE, legend="right")
10.2.3.3 Functional differences
WC from acclimation to FMT1
element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Hot_control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>10.3 What is the effect of FMT and temperature on the GM after 2 weeks caompared to acclimation?
10.3.1 CI vs CC
10.3.1.1 Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(type %in% c("Control","Treatment") & time_point %in% c("FMT2","Acclimation")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control") Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.0875) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
342.9 352.2 -165.5 330.9 29
Scaled residuals:
Min 1Q Median 3Q Max
-1.7595 -0.2361 0.2110 0.4983 1.1980
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 6.519e-12 2.553e-06
Number of obs: 35, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9512 0.1816 21.756 <2e-16 ***
typeTreatment -0.1372 0.2502 -0.548 0.583
time_pointFMT2 0.1206 0.2491 0.484 0.628
typeTreatment:time_pointFMT2 0.4149 0.3469 1.196 0.232
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typTrt t_FMT2
typeTretmnt -0.726
tim_pntFMT2 -0.729 0.529
typTr:_FMT2 0.524 -0.721 -0.718
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Modelq1n <- lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Modelq1n)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
253.6069 262.2108 -120.8034
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.268167 9.858563
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 21.455542 3.670059 16 5.846103 0.0000
typeTreatment -4.056558 5.045308 16 -0.804026 0.4332
time_pointFMT2 2.262002 4.804331 15 0.470826 0.6445
typeTreatment:time_pointFMT2 20.220508 6.684284 15 3.025082 0.0085
Correlation:
(Intr) typTrt t_FMT2
typeTreatment -0.727
time_pointFMT2 -0.697 0.507
typeTreatment:time_pointFMT2 0.501 -0.684 -0.719
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.7676771 -0.6504757 0.2922823 0.6462932 2.0120364
Number of Observations: 35
Number of Groups: 18 Phylogenetic
Model_phylo <- lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
133.508 142.1119 -60.75399
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.2974732 1.463985
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 5.233354 0.5281280 16 9.909252 0.0000
typeTreatment 0.306587 0.7258724 16 0.422370 0.6784
time_pointFMT2 -0.458617 0.7121982 15 -0.643946 0.5293
typeTreatment:time_pointFMT2 0.033201 0.9917181 15 0.033479 0.9737
Correlation:
(Intr) typTrt t_FMT2
typeTreatment -0.728
time_pointFMT2 -0.715 0.521
typeTreatment:time_pointFMT2 0.514 -0.705 -0.718
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.403145672 -0.437916442 -0.005178256 0.313089271 1.956678742
Number of Observations: 35
Number of Groups: 18 10.3.1.2 Beta diversity
Number of samples used
samples_to_keep_post5 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control") %>%
select(Tube_code) %>%
pull()
subset_meta_post5 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
subset_meta_post5$time_point<-as.factor(subset_meta_post5$time_point)
subset_meta_post5$type<-as.factor(subset_meta_post5$type)
length(samples_to_keep_post5)[1] 35Richness
richness_post5 <- as.matrix(beta_q0n$S)
richness_post5 <- as.dist(richness_post5[rownames(richness_post5) %in% samples_to_keep_post5,
colnames(richness_post5) %in% samples_to_keep_post5])
betadisper(richness_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.004294 0.0042942 0.4761 999 0.507
Residuals 33 0.297625 0.0090189
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.504
Treatment 0.49501 adonis2(richness_post5 ~ type*time_point,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 0.5174241 | 0.04775345 | 1.797587 | 0.001 |
| time_point | 1 | 0.9068994 | 0.08369841 | 3.150666 | 0.001 |
| type:time_point | 1 | 0.4878440 | 0.04502349 | 1.694822 | 0.028 |
| Residual | 31 | 8.9231562 | 0.82352465 | NA | NA |
| Total | 34 | 10.8353238 | 1.00000000 | NA | NA |
Neutral
neutral_post5 <- as.matrix(beta_q1n$S)
neutral_post5 <- as.dist(neutral_post5[rownames(neutral_post5) %in% samples_to_keep_post5,
colnames(neutral_post5) %in% samples_to_keep_post5])
betadisper(neutral_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.034255 0.034255 3.9748 999 0.049 *
Residuals 33 0.284399 0.008618
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.061
Treatment 0.054506 adonis2(neutral_post5 ~ type*time_point,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 0.6424540 | 0.06322634 | 2.589936 | 0.001 |
| time_point | 1 | 1.1894033 | 0.11705370 | 4.794863 | 0.001 |
| type:time_point | 1 | 0.6395258 | 0.06293816 | 2.578132 | 0.015 |
| Residual | 31 | 7.6897933 | 0.75678179 | NA | NA |
| Total | 34 | 10.1611764 | 1.00000000 | NA | NA |
Phylogenetic
phylo_post5 <- as.matrix(beta_q1p$S)
phylo_post5 <- as.dist(phylo_post5[rownames(phylo_post5) %in% samples_to_keep_post5,
colnames(phylo_post5) %in% samples_to_keep_post5])
betadisper(phylo_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.02670 0.026703 1.6775 999 0.226
Residuals 33 0.52532 0.015919
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.219
Treatment 0.20425 adonis2(phylo_post5 ~ type*time_point,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 0.03644585 | 0.01867131 | 0.7636208 | 0.004 |
| time_point | 1 | 0.35441680 | 0.18156873 | 7.4258113 | 0.002 |
| type:time_point | 1 | 0.08154954 | 0.04177806 | 1.7086422 | 0.175 |
| Residual | 31 | 1.47955831 | 0.75798190 | NA | NA |
| Total | 34 | 1.95197051 | 1.00000000 | NA | NA |
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post5, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post5$time_pointFMT2" = "FMT2",
"subset_meta_post5$typeTreatment" = "Cold-intervention",
"subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")
beta_richness_nmds_post5 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post5, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post5$time_pointFMT2" = "FMT2",
"subset_meta_post5$typeTreatment" = "Cold-intervention",
"subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")
beta_neutral_nmds_post5 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post5, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post5$time_pointFMT2" = "FMT2",
"subset_meta_post5$typeTreatment" = "Cold-intervention",
"subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")
beta_phylogenetic_nmds_post5 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_nmds_post5, beta_neutral_nmds_post5, beta_phylogenetic_nmds_post5, ncol=3, nrow=1, common.legend = TRUE, legend="right")
10.3.1.3 Functional differences
CC from acclimation to FMT2
element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>CI from acclimation to FMT2
significant_element<-element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Treatment") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))element_gift_sig <- element_gift %>%
select(Tube_code, all_of(intersect(
significant_element$trait,
colnames(element_gift)
))) %>%
left_join(., sample_metadata[c(1, 7, 10)], by = join_by(Tube_code == Tube_code)) %>%
filter(time_point %in% c("FMT2","Acclimation"))%>%
filter(type=="Treatment")
difference_table <- element_gift_sig %>%
select(-Tube_code, -type) %>%
group_by(time_point) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-FMT2)%>%
mutate(group_color = ifelse(Difference <0,"FMT2", "Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c('#008080',"#76b183")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
legend.title = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Function") +
ylab("Mean difference")
10.3.2 CI vs WC
10.3.2.1 Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=8))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control") Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(10.176) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
339.6 349.1 -163.8 327.6 30
Scaled residuals:
Min 1Q Median 3Q Max
-2.3741 -0.4004 0.1336 0.5764 1.4414
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 2.233e-12 1.494e-06
Number of obs: 36, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.4697 0.1104 40.482 < 2e-16 ***
typeTreatment -0.6557 0.1599 -4.101 4.11e-05 ***
time_pointFMT2 -0.2292 0.1572 -1.458 0.144814
typeTreatment:time_pointFMT2 0.7647 0.2246 3.405 0.000661 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typTrt t_FMT2
typeTretmnt -0.691
tim_pntFMT2 -0.702 0.485
typTr:_FMT2 0.492 -0.712 -0.700
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Modelq1n <- lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Modelq1n)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
255.5098 264.3042 -121.7549
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 1.317597 9.38306
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 44.60136 3.158373 16 14.121626 0.0000
typeTreatment -27.20238 4.466614 16 -6.090157 0.0000
time_pointFMT2 -12.40101 4.423217 16 -2.803618 0.0127
typeTreatment:time_pointFMT2 34.88352 6.255373 16 5.576569 0.0000
Correlation:
(Intr) typTrt t_FMT2
typeTreatment -0.707
time_pointFMT2 -0.700 0.495
typeTreatment:time_pointFMT2 0.495 -0.700 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.1054551 -0.5575213 0.2345834 0.5912103 2.0213400
Number of Observations: 36
Number of Groups: 18 Phylogenetic
Model_phylo <- lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
132.8692 141.6636 -60.43459
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7604222 1.204042
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 6.522583 0.4746882 16 13.740774 0.0000
typeTreatment -0.982643 0.6713105 16 -1.463768 0.1626
time_pointFMT2 -1.066614 0.5675911 16 -1.879195 0.0786
typeTreatment:time_pointFMT2 0.641199 0.8026950 16 0.798807 0.4361
Correlation:
(Intr) typTrt t_FMT2
typeTreatment -0.707
time_pointFMT2 -0.598 0.423
typeTreatment:time_pointFMT2 0.423 -0.598 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.90560105 -0.54341938 -0.06251794 0.39498845 1.94321495
Number of Observations: 36
Number of Groups: 18 10.3.2.2 Beta diversity
Number of samples used
samples_to_keep_post6 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
select(Tube_code) %>%
pull()
subset_meta_post6 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control")
subset_meta_post6$time_point<-as.factor(subset_meta_post6$time_point)
subset_meta_post6$type<-as.factor(subset_meta_post6$type)
length(samples_to_keep_post6)[1] 36Richness
richness_post6 <- as.matrix(beta_q0n$S)
richness_post6 <- as.dist(richness_post6[rownames(richness_post6) %in% samples_to_keep_post6,
colnames(richness_post6) %in% samples_to_keep_post6])
betadisper(richness_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.054351 0.054351 9.1422 999 0.006 **
Residuals 34 0.202132 0.005945
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.005
Treatment 0.0047271 adonis2(richness_post6 ~ type*time_point,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 1.2570893 | 0.11659822 | 4.856166 | 0.001 |
| time_point | 1 | 0.6574224 | 0.06097760 | 2.539639 | 0.001 |
| type:time_point | 1 | 0.5831993 | 0.05409322 | 2.252913 | 0.002 |
| Residual | 32 | 8.2836650 | 0.76833097 | NA | NA |
| Total | 35 | 10.7813760 | 1.00000000 | NA | NA |
Neutral
neutral_post6 <- as.matrix(beta_q1n$S)
neutral_post6 <- as.dist(neutral_post6[rownames(neutral_post6) %in% samples_to_keep_post6,
colnames(neutral_post6) %in% samples_to_keep_post6])
betadisper(neutral_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.078328 0.078328 12.006 999 0.003 **
Residuals 34 0.221814 0.006524
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.003
Treatment 0.0014541 adonis2(neutral_post6 ~ type*time_point,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 1.3255629 | 0.12619423 | 5.685213 | 0.001 |
| time_point | 1 | 0.9072453 | 0.08637020 | 3.891089 | 0.002 |
| type:time_point | 1 | 0.8102276 | 0.07713406 | 3.474989 | 0.003 |
| Residual | 32 | 7.4611123 | 0.71030151 | NA | NA |
| Total | 35 | 10.5041482 | 1.00000000 | NA | NA |
Phylogenetic
phylo_post6 <- as.matrix(beta_q1p$S)
phylo_post6 <- as.dist(phylo_post6[rownames(phylo_post6) %in% samples_to_keep_post6,
colnames(phylo_post6) %in% samples_to_keep_post6])
betadisper(phylo_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04938 0.049381 3.8729 999 0.036 *
Residuals 34 0.43351 0.012750
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.036
Treatment 0.057271 adonis2(phylo_post6 ~ type*time_point,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 1 | 0.1600877 | 0.08395907 | 3.885816 | 0.001 |
| time_point | 1 | 0.2882491 | 0.15117420 | 6.996685 | 0.001 |
| type:time_point | 1 | 0.1400636 | 0.07345727 | 3.399769 | 0.030 |
| Residual | 32 | 1.3183347 | 0.69140946 | NA | NA |
| Total | 35 | 1.9067351 | 1.00000000 | NA | NA |
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post6, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post6$time_pointFMT2" = "FMT2",
"subset_meta_post6$typeTreatment" = "Cold-intervention",
"subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")
beta_richness_nmds_post6 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post6, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post6$time_pointFMT2" = "FMT2",
"subset_meta_post6$typeTreatment" = "Cold-intervention",
"subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")
beta_neutral_nmds_post6 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post6, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post6$time_pointFMT2" = "FMT2",
"subset_meta_post6$typeTreatment" = "Cold-intervention",
"subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")
beta_phylogenetic_nmds_post6 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_nmds_post6, beta_neutral_nmds_post6, beta_phylogenetic_nmds_post6, ncol=3, nrow=1, common.legend = TRUE, legend="right")
10.3.2.3 Functional differences
WC from acclimation to FMT2
element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Hot_control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>